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> : 

Q , Exploiting an approach similar to the R-matrix theory, the diffusion Monte Carlo method is em- 

. ployed to compute phase shifts and threshold cross sections for the elastic scattering of o-positronium 

' off light atoms. Results are obtained for Ps-H and Ps-He as representative cases of open and closed 

\^ ^ shell targets. The method allows for an exact treatment of both correlation and exchange interac- 

■ tions, and represents the most promising approach to deal with these effects in more complicated 

targets. In particular the Ps-He threshold cross section, computed in a many body framework for 
the first time, represents a standard by which past and future numerical and experimental estimates 
can be judged. 
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' Positronium (Ps) scattering off atomic and molecular targets has an overwhelming importance if an understanding 
^ ' of the interaction mechanism between an overthermal Ps and a condensed matter environment is required [1]. For 
Q instance, by means of elastic and inelastic cross sections, it may be possible to model energy transfers from Ps to the 
' surroundings or to describe the trapping of Ps in a free volume cavity. Despite its long history [2, 3, 4, 5, 6, 7], and 
^»-^, even in the case of light atoms, some quantitative aspects of the process still remain controversial and have recently 
' been addressed by a number of authors, both experimentally [8, 9, 10, 11] and theoretically [12, 13, 14, 15, 16]. 
OhI, From the computational point of view, the difficulties with which almost every method is faced are related to the 
composite nature of both target and projectile. As a consequence, sensible results can be obtained only if correlation 
^-H ' and exchange effects are properly treated. Moreover, an accurate description of the correlation effects between the 
^ : target electrons and the positron in the Ps is important in computing "pick off' annihilation rates. These effects have 
only recently been computed consistently for the case of positronium scattering off hydrogen and positronium atoms 
, [14, 17]. However, the full ah initio treatment (i.e. without the use of exchange or correlation model potentials) of 
' systems with more than two electrons still represents a formidable task. A glance to the recent literature on bound 
I j systems containing a positron reveals a practically identical situation with only small numbers of electrons treated 
I • explicitly. In this context, it has been shown by several authors [18, 19] that flexible and accurate computational 
[ techniques for small and medium size systems are provided by the family of quantum Monte Carlo (QMC) methods. 
. Among them, the diffusion Monte Carlo (DMC) scheme represents the most powerful approach to study strongly 
^ ■ correlated systems thanks to its ability to sample a distribution proportional to the exact ground state wave function 
' ^ , of a given Hamiltonian. Where fermions are concerned, the antisymmetric nature of the wave function and its 
consequent non-positiveness imply the introduction of a bias known as nodal error. For a given energy e, the nodal 
, error Ae, which disappears if the nodal surfaces of the exact wave function arc known, has a value which commonly 
spans the range Ae/e G [10 ^, 10 ^] [20]. Unless otherwise specified, the following results have been computed in the 
^ ■ DMC framework. 

The application of QMC methods to scattering problems was independently proposed in the eighties in two pi- 
oneering papers by two groups in the field of nuclear physics [21, 22]. Their ideas have been recently applied to 
the exciton-exciton scattering problem [17], thus providing the first accurate calculation for the Ps-Ps system. This 
approach, which closely resembles the original idea behind the R-matrix theory of Wigner [23], is briefly summarized 
in what follows for the case of an elastic collision. We define yab = R-A ^ R-s as the relative position of the centers 
of mass of A and B and p and jj, as their relative momentum and reduced mass. The choice of a boundary tab = Ti- 

2 

under the condition Vint{T^) ^ f^' allows one to approximate the wave function in the region tab > with the 
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asymptotically exact form 



^ = A 



(1) 



A is the antisymmetrization operator, sa and the internal coordinates of the two separate fragments, "^a and 
their wave functions and Yim and the angular and radial functions describing the dynamics of the relative motion 
of the two centers of mass. 

The stationary form of can be expressed as 

^livAB) = Mptab) + Si{p)Oi{prAB) (2) 
where X; and Oi are Hankel functions, and Si{p) is the scattering matrix. Here, p is connected to the total energy 

2 

by Etot = 2^ + Ea + Eb where the computation of Ea and Eb-, both being bound state energies, can be exactly 
performed employing the DMC method. The condition = ^ ^md a corresponding value of p (i.e. Etot) are 

enough to precisely state the value of Si{p) in Eq.(2). The link between B and p is contained in the dynamics of the 
interior region and therefore the computation of Eioi in this region, under the same boundary condition on tab = 
(i.e. imposing the continuity of the logarithmic derivative), provides the value of Si. Eventually the scattering matrix 
assumes the form 

OlipTl) - BOiipU) 

Hence, as long as the boundary condition is exactly controlled, all the relevant information is contained in the interior 

region. In this subspace the wave function is normalizable and the energy, now paramctrically dependent on TZ, can 
be studied by means of one of the QMC techniques e.g. VMC, DMC or one of their variants. In order to easily fulfill 
the above boundary condition in the DMC framework, one can choose the value of B to be infinite. Prom a physical 
point of view, this choice corresponds to add a rigid wall located at TZ along the distance between the target and 
projectile centers of mass. Results will be presented in terms of phase shift 6i{p) which is defined by Si = e^**' and 
can be expressed, starting from Eq.(3) and taking the B — > oo limit, as 

where ji and n; are respectively the spherical Bcssel and Neumann fimctions. 

Before going on we feel that it is worth stressing two important points. First, one has to satisfy the condition 

2 

y <^ f-, SO an upper limit to the sphere radius TZ does not exist, while it cannot be chosen smaller than some 
unfortunately not well specified threshold value. This imposes an upper limit to the relative kinetic energy. For Ps 
scattering off neutral atoms, the interaction potential between the target and the projectile dies off as l/r\g. This 
short range potential allows the use of fairly small radii, a possibility not necessarily available for different colliding 
fragments. Secondly, since DMC samples the lowest energy state, it cannot be applied to scattering problems in 
presence of a bound state. 

Both these two points highlight the importance that studying excited states could have within this approach. The 
possibility of raising the energy whilst keeping the surface constrain fixed can, in principle, allow the study of every 
system at any energy. 

In this work, we applied the presented technique to the S-wave scattering of positronium off hydrogen and helium 
dynamically described by the full Hamiltonian 

N N N 

^ i=i ^ i=i i>j ^'-j i 

where i and j refer to electrons, p to the positron and Z to the nuclear charge of the atom. The spatial part of the 
scattering wave function has been chosen to be of the form 

* = 0[M/A(sA)*p.(rip)^^^^</>j(s/)] (6) 

rpsA 

where ^a, '^Ps, and <I> have the same meaning as in Eq. (1). is a Jastrow factor for all the pairs of particles 
belonging to different fragments, s/ is the set of distances for these pairs and O is the appropriate symmetry operator 



3 



TABLE I: Scattering lengths (bohr) for Ps-H scattering. 




QMC 


Previous results 


S=0 


4.36(2) 


4.3 [14] , 3.49 [13] , 5.20 [12] 


S=l 


2.24(1) 


2.2 [14] , 2.46 [13] , 2.45 [12] 



built according to Young diagrams. In the Ps-H case, the exact internal wave function of both fragments is known 
and O has the form 

= 1 + (-l)^Pi2 (7) 

where 5 (0 or 1) is the spin momentum of the state and P12 the permutation operator between the two electrons. The 
5 = space part of the wave function is everywhere positive, while the nodal surface for the S = 1 state is exactly 
provided by the action of O. Under this condition the energy can be computed by DMC without any approximation. 
The singlet state supports a bound state and, as said above, it is thus necessary to exploit an excited state technique. 
Evidently, in order to get the required scattering information, the boundary condition at the surface for every state 
must be controlled. As shown in Ref [17] the Correlation Function DMC method[24] with our choice of B accomplishes 
automatically this requirement and it will therefore be used in the following. Detailed descriptions of this method 
are out of the scope of the present work and can be extensively found in the literature [24, 25]. In this respect, we 
limit ourselves to show the convergence for one calculation in Figure 1. All the simulations for the triplet state of 
Ps-H were carried out using a time step of 0.01 hartrce"'^, 2000 walkers, and a total of 100 blocks of 10000 steps each. 
Simulations for the singlet states were performed employing 2000 configurations, a time step of 0.01 hartree""'^, and 
a grand total of 10000 decorrelated Euclidean time evolutions. Low cwigy phase shifts for both 5 = and 5 = 1 
systems are shown in Figure 2. Scattering lengths are reported in Table I, together with three independent estimates 
of the same quantities very recently calculated [12, 13, 14]. Whereas all of them agree in assigning the value for the 
triplet state, tlicax; appears to be some controversy where the singlet state is concerned. Our results, which we believe 
to be statistically exact, are very close to the values proposed in Ref [14] suggesting these could be safely considered as 
a definitive estimate. This fact can be also taken as a strong proof of the reliability of the method we are employing, 
as well as of the full-electron Stochastic Variational Method approach used in Ref [14]. 

With this premise, we now address the more debated problem of positronium scattering off helium. Before discussing 
our computed quantities for this process, it is worth noting that the experimental measurements of the threshold value 
of the cross section span almost an entire order of magnitude [8, 11]. The most recent theoretical estimates, obtained 
by different computational schemes, do not single out one of these as the correct one. The primary reason for this 
failure is the small size of the cross section and a consequent large fractional error associated with any approximation. 
To make a more direct comparison with experiments. Table II shows the most recently calculated and measured 
threshold cross sections. 

In the present study, the system is treated with a genuinely many-body technique and no physical approximation 
have been made prior to the numerical simulation. The absolute freedom one has in choosing the analytical form of 
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FIG. 1: Energy decay for the first four states of Ps-H system with S=0 and 71=15 bohr. Note that the energy of the ground 
state is correctly less than -0.75 hartree and coincident with that of the bound state of PsH, -0.78919 hartree. 



FIG. 2: Phase shift for Ps-H S-wave elastic scattering with total electron spin S=l and S=0. Momentum in atomic units and 
phase shift in radians. 
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FIG. 3: Phase shift for Ps-He S-wave elastic scattering (S=l/2). Momentum in atomic units and phase shift in radians. 



the wave function in QMC methods allows us to employ the following explicitly correlated form for '^He 

*He = expf^^i!^^ + %tM_,^l!^) (8) 

which has a DMC energy statistically exact. Moreover, the only Young diagram compatible with the choice of an 
helium atom in its ground state {S = 0) gives the following form for O: 

= (1 + Pl2)(l-Pl3) (9) 

Simulations for this system were characterized by a time step of 0.005 hartree^^, 4000 walkers, and a total of 130 
blocks of 25000 steps each. The numerical results for the phase shift are shown in Figure 3. 

The value of the computed scattering length is 1.4046(6) bohr with a corresponding threshold cross section of 
7.892(2)7raQ. Comparison with numerical estimates and experimental results curiously shows this value to lie in a 
'neutral' zone, intermediate amongst the most recently proposed values (Table II). The nodal error, being the only 



TABLE II: Scattering threshold cross section (tt bohr^) for He-Ps. 



QAIC 


Expt'riiiioutal 


I'roviou.s calculations 


7.892(2) 


13(4) [8] 


13.2 [16] 




2.6(5) [11] 


3.10 [15] 






10.4 [14] 
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approximation introduced, deserves some comments. General considerations [26] show this bias on the phase shift to 
be always negative and proportional to V^^/"^ where V is the sphere volume. As a result of this, our scattering length 
could be slightly lower than the exact one. More quantitatively, one can observe that in the interaction region the 
employed function closely resembles the functional form used in bound state calculations on similar systems, for which 
the nodal error roughly equals 10~^ hartree [19]. In the rest of the simulation volume the nodes of the trial wave 
function are practically exact thanks to the validity of Eq. 1. For such a reason we expect a bias on the energy of the 
same order of magnitude as the one in bound state calculations. If so, the nodal error would turn out to be smaller 
than the statistical fluctuations of our energy values, therefore warranting the statistical exactness of our results. 

Among the numerous applications directly derivable from this method, we would like to emphasize that the possibil- 
ity to sample the exact particle distributions in configurational space may allow one to obtain an effective interaction 
potential between Ps and a given atom or molecule. This potential, where all the physical effects are correctly ac- 
counted for, could be successively used to simulate Ps in condensed phases as molecular crystals and liquids. Moreover, 
it could also help in defining the preferential spatial location where the Ps positron would annihilate during a "pick 
off' annihilation event, so that the interplay between the theoretical and the experimental results may enhance the 
diagnostic role played by Ps in condensed matter science. The DMC method is also suitable, as demonstrated in a 
slightly different context [27], for the computation of Zeff, and its formal extension to reactive processes (inelastic 
are still excluded) has been known since the seminal work of Alhassid and Koonin [21]. 

The authors are in debt to Dr. Jim Mitroy for many helpful comments and discussions on methodological issues 
and positronium physics. Financial support from the Universita of Milano is also acknowledged. 
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